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Charge order has emerged as a generic feature of doped cuprates, leading to important questions 
about its origin and its relation to superconductivity. Recent experiments on two classes of hole 
doped cuprates indicate a novel d-wave symmetry for the order. These were motivated by earlier 
spin fluctuation theoretical studies based on an expansion about hot spots in the Brillouin zone that 
indicated such order would be competitive with d-wave superconductivity. Here, we reexamine this 
problem by solving strong coupling equations in the full Brillouin zone for experimentally relevant 
parameters. We find that bond-oriented order, as seen experimentally, is strongly suppressed. We 
also include coupling to Big phonons and do not see any qualitative change. Our results argue 
against an itinerant model for the charge order, implying instead that such order is likely due to 
Coulombic phase separation of the doped holes. 
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Following earlier NMR studies^, recent x-ray experi¬ 
ments on underdoped cuprates have detected short range 
charge density wave (CDW) order with a period of 3 to 
5 lattice spacings in YBCO^r— , and in Bi based^i^ and 
Hg based cuprates^. For YBCO, the doping range where 
the charge order has been observed^ coincides with the 
doping range where quantum oscillation experiments de¬ 
tect reconstruction of the Fermi surface^. In conven¬ 
tional CDW systems, the charge order is thought to 
have s-wave symmetryiS. In contrast, scanning tunnel¬ 
ing microscopyii and resonant soft x-ray scatterin g^^d^ 
data have revealed a novel d-wave symmetry, where the 
two oxygen ions in a Cu02 unit are out of phase. This 
charge order differs from the more robust stripe order 
seen earlier in La based compoundsiA2S, which appears 
to have s-wave symmetry instead^^. In all cases, though, 
the wavevector is oriented along the bond directioiiii. 

The search for d-wave symmetry was motivated by 
earlier theoretical studies of Metlitski and Sachde v^^d^ . 
They have shown that charge order is competitive with 
d-wave superconductivity in a spin fluctuation model. 
This instability has a d-wave form factor, with a diagonal 
wavevector that spans Fermi surface points (hot spots) 
that intersect the antiferromagnetic zone boundary of the 
undoped phase (Fig. [T|). The subsequent observation of 
charge order in YBCO motivated a number of follow up 
studies^S^^. Most of these studies are either based on 
an expansion around the hot spots, with the Fermi sur¬ 
face curvature treated as a perturbatio n^°d^d^ , or rely on 
a weak coupling approximation. So, the question arises 
whether these results survive in a strong coupling treat¬ 
ment where these approximations are not made. 

Here, we solve the strong coupling instability equa¬ 
tion for the charge order in the entire Brillouin zone in¬ 
cluding the full momentum and frequency dependence 
of the bosonic and fermionic spectra. This formalism 
has been used in the past to study d-wave supercon¬ 
ductivity originating from spin fluctuations^^i^. It has 
also been used to study instabilities in the particle-hole 


channel^l. Recently, this formalism was used by us to 
examine the effect of the pseudogap on spin-fluctuation 
mediated pairing^S. The one approximation we make is 
that the bosonic and fermionic spectra are taken from 
experiment rather than self-consistently calculated. We 
find that bond-centered charge order is completely sup¬ 
pressed, and inclusion of the dressed fermion Green’s 
function additionally suppresses diagonal charge order, 
with the only robust order in this model for experimen¬ 
tally relevant parameters being d-wave superconductiv¬ 
ity. 

Our starting point is the linearized equation for the 
anomalous self energy in the particle-hole channel 



FIG. 1: (Color online) Fermi surface for the tight binding 
dispersion considered in this work. The dashed lines show 
the antiferromagnetic Brillouin zone and the dotted lines the 
structural one. Filled circles denote the hot spots and filled 
squares the antinodal points. The wavevectors studied here 
are indicated by the arrows. 
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T ^ V{k-k',iuJn-iuJm)'P'^{k',iuJm)^^{k',iuJm) 

k' ,UJrn 

= A$'3(fc,^c^„). (1) 


Here Q is the ordering vector, V is the interaction and 

is the CDW particle-hole kernel 

V^{k',iu^m) = G{k' - ^^iuJm)G{k' + (2) 

where G is the fermion Green’s function. Because of the 
complexity of the strong coupling equations, as a first 
approximation, we do not include the frequency depen¬ 
dence of the fermion self-energy, and thus set G(k, iuj) = 
(iw —For ^fc, we consider a renormalied dispersion 
that fits low energy ARPES data for Bi 2 Sr 2 CaCu 2 08+5 
(tb2 dispersion of Ref. [l|). Later in the paper, we will 
include the fermion self-energy as well. The interaction 
assumed here is 




dx V"{k, x) 

TT — X 


(3) 


where V" is proportional to the imaginary part of the dy¬ 
namic spin susceptibility. We consider the phenomeno¬ 
logical formas 


V"{k,n) 

Xk 


2 

l^dsfXQ 


Xk^lf + 12 ^ 


{^af/o) ^ + 2 + cos k^a + cos kytt 


(4) 


where gsf is the spin-fermion coupling constant, ^af is 
the antiferromagnetic correlation length, XQ is the static 
susceptibility at Qaf = (7r/a,7r/a) and 12^/ is the char¬ 
acteristic spin-fluctuation energy scale. Because of the 
1/12 decay of V", we impose a frequency cut-off, 12^. We 
use 12c =300 meV, 12^/ =100 meV, gIfXQ =0.9 eV and 
^af = 2a, where a is the lattice constant. The values of 
^AF and 12s/ are motivated from inelastic neutron scat¬ 
tering studies of the magnetic excitations of underdoped 
YBCO near ■ 12c is motivated from recent res¬ 

onant inelastic x-ray scattering studies of the higher en¬ 
ergy excitations away from Qap ^^'^'^ - The value of g'^fXQ 
was chosen to obtain a superconducting Tc of 50 K, typi¬ 
cal for underdoped YBCO. At the transition, the leading 
eigenvalue A in Eq. m reaches 1, with its eigenvector 
giving the structure of the CDW order parameter. 

We consider the CDW instabilities for the diagonal 
CDW case (CDW-diag) with ordering vector {Qhs, Qhs), 
and for the bond oriented case (CDW-x) with vectors 
(Qan,0) and {Qhs,0), as shown in Fig. 1. We perform 
our calculations with a 0.027r/a momentum grid with 32 
Matsubara frequencies, which is sufficient for convergence 
of the eigenvalues for the temperature range studied here 
(see supplementary information). Fig. [5] shows the tem¬ 
perature dependence of the leading eigenvalue for the dif¬ 
ferent CDW cases along with the d-wave superconducting 
case. 

As expected, the eigenvalue for the superconducting 
case exhibits a logarithmic divergence in T. This is 



FIG. 2: (Color online) Temperature dependence of the lead¬ 
ing eigenvalues for the superconducting d-wave state (SC), 
the diagonal CDW state (CDW-diag) and the bond oriented 
CDW state (CDW-x), in an approximation where G is based 
on a renormalized dispersion taken from ARPES data. The 
ordering vector for CDW-diag is {Qhs, Qhs)- For CDW-x, two 
vectors were considered: (Qan,0) and (Qhs,0). 


present as well for CDW-diag order, though we find it 
to be significantly reduced relative to the superconduct¬ 
ing one. The CDW-diag state has d-wave (Big) symme¬ 
try with a momentum dependence that is well described 
by cos(fca;a) — cos(kya) as can be seen in Fig. [3] (a). In¬ 
creasing the antiferromagnetic coherence length doesn’t 
change our findings. The CDW-diag instability becomes 
stronger with longer ^af , but it always remains subdom¬ 
inant relative to d-wave superconductivity (see supple¬ 
mentary material). 

We now focus on the bond-oriented CDW states, since 
there is no experimental evidence for diagonal-oriented 
order. The T dependence of the leading eigenvalues for 
CDW-x are also plotted in Fig. [21 They are almost iden¬ 
tical for vectors (Q/is,0) and (QanjO)- The eigenvalues 
vary little with temperature, showing no evidence for a 
log in the temperature range studied. This is one of our 
main results, and differs from an analytic approximation 
quoted in earlier worb^i. The lack of a log divergence 
in our case is likely due to the curvature coming from 
the experimentally determined Fermi surface. Moreover, 
our bosonic spectrum is based on inelastic neutron and 
x-ray scattering measurements. Hence, we are not in the 
extreme quantum critical regime considered in Ref. 

We next look at the structure of the CDW-x state. 
Fig. 121(b) shows the momentum dependence of the eigen¬ 
vector corresponding to the leading eigenvalue, and can 
be well fit by a sum of a constant (s), cos{kxa) + cos{kya) 
(s') and cos{kxa) — cos(fcga) (d), with the d-wave compo¬ 
nent dominant, consistent with earlier studies^!. 

We now turn to the influence of the fermion self-energy. 
The fermion self-energy in a lowest order approximation 
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FIG. 3: (Color online) The momentum dependence of the eigenvectors at the lowest Matsubara frequency (r=40 K) corre¬ 
sponding to leading eigenvalues (Fig. 2) for the (a) CDW-diag and (b) CDW-x states. For CDW-x, the ordering vector is 
(Qan,0), though the eigenvector for (Qhs,0) is similar. A normalization condition of j, = 1 is employed. 


2 . 


1 . 


0 . 


FIG. 4: (Color online) Temperature dependence of the leading 
eigenvalues for the superconducting d-wave state (SC), the di¬ 
agonal CDW state (CDW-diag) and the bond oriented CDW 
state (CDW-x), obtained using dressed Green’s functions. 
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I](fc,fw) = TV{k — k',iu} — iu}')Go{k' (5) 

k' ,uj' 

= V{k-k',iu;-iu;') ( 6 ) 

k' ,uj' ^ 

where Go is the bare Green’s function. Since our previous 
was based on ARPES data, we scale it by a factor of 
two to get an approximation to the bare dispersion, with 
a bare nodal Fermi velocity =3.2eV-A. For the spin- 
fluctuation interaction, we use g^fXQ = 3.2eV and keep 
the rest of the parameters the same as before. This cou¬ 
pling strength renormalizes the nodal Fermi velocity to 
1.8eV-A, which is comparable to the experimental value. 


To keep the shape of the Fermi surface intact, we drop 
the term proportional to ^k' in the self-energy calcula¬ 
tion (Eq. (5]). We now use this dressed Green’s function 
{G~^{k,ii-o) = iui — ^k — E(fc,fw)) for the instability anal¬ 
ysis in Eq.[T] Its effect (Fig. Id]) is to additionally suppress 
the diagonal charge order (presumably due to the energy 
smearing of G). As a consequence, for experimentally 
relevant parameters, only d-wave superconductivity re¬ 
mains as a robust instability in this model. 

Our results cast doubt on an itinerant spin fluctuation 
mediated origin for the observed charge order. On the 
other hand, the dependence of the observed wavevector 
on doping^ is suggestive that the Fermi surface is playing 
some role as in classic CDW systems. This is in contrast 
to the La based cuprates whose doping dependence is op¬ 
posite to this, as would be expected from a real space pic¬ 
ture where the wavevector is proportional to the doping. 
In classic CDW systems like ^R-NbSe 2 , phonons play an 
important role^^. Hence, one could think that this might 
be the case in the cuprates as well, where anomalies have 
been seen in both optio^ and acoustic^ phonon modes 
near the charge ordering vector. It is interesting to note 
that Big phonon modes have been postulated to be re¬ 
sponsible for dispersion anomalies seen in photoemission 
near the antinodes^, and perhaps their d-wave symmetry 
is related to that of the charge order. In support of this, 
several theoretical studies have suggested that coupling 
of the electrons to such modes can cause d-wave charge 
orde r^^4° . We examine this idea by including the phonon 
mediated interaction along with the spin-fluctuation me¬ 
diated interaction in the CDW instability equation. Un¬ 
der the assumption that the electron-phonon interaction 
is modest, we do not find that its addition has sufficient 
strength to cause a CDW instability either (see supple¬ 
mentary material). 

None of the above means that spin mediated interac¬ 
tions are irrelevant for the charge order, but our results 
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indicate that an itinerant treatment of the problem in a 
parameter range appropriate to experiment is not ade¬ 
quate to capture the resulting physics. In that context, 
we point to real space treatments of the t — J model like 
density matrix renormalization group (DMRG) which in¬ 
dicate a strong tendency to bond-centered charge order 
in the underdoped regime^, with the d-wave structure 
likely due to Coulomb repulsion which acts to enforce 
charge neutrality in each unit cell. 

Note added in proof. After the original submission of 
this work, a critique of it has been offered by Wang and 
Chubuko v^^d^ . Their primary criticism revolves around 
two points - that in the original submission, we neglected 
the fermion self-energy that they argue lessens the detri¬ 
mental impact of the Fermi surface curvature, and that 
our value of flsf is too small by an order of magnitude. In 
regards to the first point, using dressed Green’s functions, 
we now hnd strong suppression of the charge order, even 
for diagonal wavevectors. In regards to the second point, 
our value of Qsf is based on experiment, and an order of 


magnitude larger value would lead to a spin fluctuation 
energy scale at the F point of the zone of 4 eV, in gross 
disagreement with inelastic x-ray scattering data. We 
also comment that the logarithmic temperature diver¬ 
gence of the eigenvalue derived for bond centered charge 
order in their work is based on analytic approximations 
which we feel are not valid for experimentally realistic 
parameters. Gertainly, we find no evidence for a log in 
our own numerical studies, even for a ^af as large as lOo 
(see supplementary material). As an additional note, our 
results are consistent with recent findings based on a real- 
space version of Eliashberg theory^. 
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Supplementary material 

Effect of a larger antiferromagnetic coherence length 



Supplementary Fig. 1: Temperature depeirdence of the lead¬ 
ing eigenvalues for the superconducting (SC), diagonal CDW 
(CDW-diag) and bond-oriented CDW (CDW-x) states, for 
^AF = 5a and 10a, in an approximation where G is based on 
a renormalized dispersion taken from ARPES data. 

Supplementary Fig. [T] shows the effect of larger anti¬ 
ferromagnetic correlation lengths. We use the same pa¬ 
rameters for interactions as for ^af = 2a, except for the 
overall prefactor g^fXQ^ which we set to 0.49 eV and 0.35 
eV for ^af = 5a and 10a, respectively, in order to obtain 
a superconducting Tc of 50 K. As can be seen, increasing 
the antiferromagnetic correlation length leads to similar 
results to those presented in the main text for ^af = 2a 
(Fig. 2). 

Self-energy effects 

Now, we include the effect of the fermion self-energy 
on the temperature dependence of the eigenvalues. Sup¬ 


plementary Fig. [2] shows the temperature dependence of 
the eigenvalues for larger antiferromagnetic correlation 
lengths. For this study we use g^fXQ = 3.2eV and keep 



T/K T/K 


Supplementary Fig. 2: Temperature dependence of the lead¬ 
ing eigenvalues for the superconducting (SC), diagonal CDW 
(CDW-diag) and bond-oriented CDW (CDW-x) states, for 
^AF = 5a and 10a, obtained using dressed Green’s functions. 
The left column shows the temperature dependence of the 
eigenvalues when the Fermi surface is kept unchanged in the 
full Green’s function dropped in the numerator of Eq. 6), 
whereas for the figures in the right column, the full expression 
for G in Eq. 6 is used. 

all other parameters the same as in the main text. 


Effect of the band structure on the CDW order 

We study a different fermionic dispersion where Qan 
and Qhs are further separated than in the previous 
dispersion to test how this influences the solutions at 
(Qan,0) and (QhsjO). Supplementary Fig. [3] shows the 
Fermi surface for this alternate band structure and the 
resulting temperature dependence of the leading eigen¬ 
values for various CDW states. Here we use ^af = 5a 
and keep the rest of the parameters for the interaction as 
before, except again we adjust g^fXQ to 0.73 eV in order 
to obtain a superconducting Tc of 50 K. The fermionic 
dispersion is 

^{kx,ky) = 0 A {cos kxtt + cos kya) — 0.32 cos k^a cos kyU 
-1-0.04 (cos 2fca;a-I-cos 2A:ya) — 0.15 (7) 

where all energy scales are in eV. For this dispersion, 
we can clearly see that the eigenvalues of the CDW-x 
state for Qhs and Qan are quite different, though neither 
exhibit a log. We have also tested the tbl and tb4 disper¬ 
sions of Ref. [l| and do not find any qualitative difference 
in our conclusions. 
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Supplementary Fig. 3: Fermi surface for the altenrate 
fermionic dispersion is shown in the upper panel. The lower 
one shows the temperature dependence of the leadiirg eigen¬ 
value for the superconducting (SC), diagonal CDW (CDW- 
diag), and bond-oriented CDW (CDW-x) states for ordering 
vectors Qha and Qan, in an approximation where G is based 
on a renormalized dispersion taken from ARPES data. 


Effect of the Big phonon mode on the CDW order 


We consider the Big phonon mode which involves an 
out of phase vibration of adjacent oxygen atoms along 
the c-axis. In general, the electron-phonon matrix ele¬ 
ment is sensitive to the details of the underlying elec¬ 
tronic structure. We use the same fermion dispersion as 
in the main text, and for the electron-phonon matrix el¬ 
ement we consider that based on a three-orbital model 
for the cuprates^. The effective interaction reads (with 
Vs f the spin-fluctuation mediated interaction of the main 
text) 




eff 


= Vsf-g{k'-Q12, k 
2flBlg 


Ql2)g*{k' + Ql2,k + Ql2) 
( 8 ) 



Supplementary Fig. 4: The temperature dependence of the 
CDW eigenvalues with both spin-fluctuation and phonon me¬ 
diated interactions, in an approximation where G is based 
on a renormalized dispersion taken from ARPES data. The 
dashed line is the eigenvalue for the superconducting state, 
which gets slightly enhanced with the inclusion of the phonon 
mediated interaction. 


where ilsig is the Big phonon mode energy which we 
take to be 341 cm“^, and g{k, k') is the electron-phonon 
matrix element which for this mode has a strong depen¬ 
dence on k and k—. For the overall magnitude of g con¬ 
sidered here, the Fermi velocity of 3.2 eV-A^ along the 
zone diagonal from our tight binding dispersion would 
be renormalized to 2.7 eV-A. Results are show in Sup¬ 
plementary Fig. 4. 


Convergence of eigenvalues 

We have performed our calculations on a 101x101 mo¬ 
mentum grid, which corresponds to a 0.027r/a momentum 
resolution. We checked different grid sizes and find that 
this grid is sufficient for convergence of the eigenvalues. 
Supplementary Fig. [5] shows the variation of eigenvalues 
as a function of the grid size. We also find that the num¬ 
ber of Matsubara frequencies used in our calculation is 
sufficient for convergence of the eigenvalues in the tem¬ 
perature range that we study. 


Eigenvector projection on the Fermi surface 


In Supplementary Fig. 6 we show the projection of the 
eigenvector on the Fermi surface for the CDW-diag and 
CDW-x states. The full Brillouin zone dependence of 
these eigenvectors are shown in the main text (Fig. 3). 
On the Fermi surface, the CDW-diag state is well de¬ 
scribed by a cos2(/) function, where (j) is the angle along 
the Fermi surface. The bond oriented CDW-x state can 
be fit with a + b cos 2(j) + c cos A(j), where these terms rep- 





























resent the s, d, and s' components, respectively, with the 
d component the dominant contribution. 



Supplementary Fig. 5: The dependence of the maximum 
eigenvalue (as in Fig. 2) on the grid size Nk x Nk- These 
eigenvalues correspond to the CDW-diag state at 40 K. In 
each case, we use 16 Matsubara frequencies. 
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Supplementary Fig. 6: The dependence of the eigenvector cor¬ 
responding to the leading eigenvalue (Fig. 2) for the ordering 
vectors (Qhs,Qhs), and (Qan,0), projected on the Fermi sur¬ 
face, where tj) =0 corresponds to the antinodal point. This is 
shown for the lowest Matsubara frequency at T=40 K. 
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